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ABSTRACT 


A low Reynolds number k - t turbulence model has been presented which yields ac- 
curate predictions of the kinetic energy near the wall. The model is validated with the 
experimental channel flow data of Kreplin and Eckelmann. The predictions are also com- 
pared with earlier results from direct simulation of turbulent channel flow. The model is 
especially useful for internal flows where the inflow boundary condition on e is not generally 
easily prescribed. The model partly derives from some observations based on earlier direct 
simulation results of near-wall turbulence. Very near the wall, limiting behavior on the 
turbulence kinetic energy and the length scale is imposed. This minimizes the sensitivity of 
the model to the inflow length scale and therefore the inflow dissipation rate prescription. 
The dissipation rate equation is modified near the wall so that the molecular diffusion and 
the destruction terms balance each other uniformly as the wall is approached. A new / M 
function is prescribed that is consistent with the limiting behavior of k and e near the wall. 
Alternate boundary conditions on e have been derived. 

The low Reynolds number turbulence model together with an existing curvature cor- 
rection appropriate to spinning cylinder flows has been used to simulate the flow in a 
U-bend with the same radius of curvature as the Space Shuttle Main Engine (SSME) 
Turn-Around Duct (TAD). The present computations indicate a space varying curvature 
correction parameter as opposed to a constant parameter as used in the spinning cylinder 
flows. The predictions show high levels of turbulence on the concave wall as expected. 
Comparison with limited available experimental data is made. The comparison is favor- 
able, but detailed experimental data is needed to further improve the curvature model. 


INTRODUCTION 

The numerical solution of turbulent flow in internal flow systems is becoming in- 
creasingly necessary both as a pre-design and a post-design engineering tool. Amongst 
numerous flow applications, the internal flow in a duct with a sharp curvature has been 
one of the most difficult to handle. One of the main reasons for this difficulty is a paucity 
of the appropriate experimental data and the lack of appropriate turbulence models for 
flows with very strong curvature effects; the two not being necessarily mutually exclusive. 

Achieving a solution of turbulent flow in a 180° bend is complicated by the fact that 
the strain rate field undergoes drastic changes as the flow turns around the bend; faster 
the rate of curvature, more severe the change in the strain rates. This curvature-induced 
effect can be taken into account in a turbulence model by, e.g., incorporating an appro- 
priate curvature correction term into a k - t turbulence model. Various such corrections 
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have been proposed [1,2]. Most of these corrections owe their origin to Bradshaw’s analogy 
[3,4] of such flows to the buoyant flows and the definition of a dimensional group called 
the “rotational” or “curved flow” Richardson number. Mixing length is defined to be a 
linear function of this rotational Richardson number and the correction term is based on 
this number. Some authors have incorporated the Richardson number based correction 
term into the turbulence kinetic energy, k, equation and others into the dissipation rate, e, 
equation. In one such model 1 , the authors define a term proportional to the angular mo- 
mentum so that for a positive rate of change of angular momentum in the radial direction, 
the Richardson number remains positive and vice-versa. They include this term in the 
e equation such that the dissipation rate increases near the convex wall and it decreases 
as the concave wall is approached. This results in an increase in kinetic energy near the 
concave wall and a decrease in k near the convex wall as compared with a straight channel 
case. This correction term involves a model parameter which appears as a coefficient and 
its magnitude has been validated based on the flow over spinning cylinders 1 . 

In this study, the curvature correction model 1 has been modified to allow for an 
angular variation of this constant. This is important because curvature effects start to 
build up as the flow turns around the bend and are not representative of the situation that 
exists in flow over spinning cylinders. A form of this angular variation is presented in this 
study. In addition, cross-flow variation of this constant from convex to concave walls is 
also indicated. 

The k - e turbulence solver 5 , KEM code, is used in the present study. The curvature 
correction to the e equation is discussed below. The KEM code is cast in generalized 
coordinates in three dimensions and includes the diffusion terms in all the three directions. 
It can be used in conjunction with a compressible or an incompressible flow solver. In the 
present study, it is used in conjunction with an incompressible flow solver, INS3D code 6 , 
which has evolved from the two-dimensional 7 and axi-symmetric 8 flow solvers respectively, 
and is based on the concept of artificial compressibility 9 . Both KEM and INS3D codes 
use the approximate factorization algorithm 10 . The former solves a 2x2 system implicitly, 
while the latter a 4x4 system. These two codes were also used in the computation of 
a large separated turbulent flow region in a channel with a sudden expansion 11 and in 
the simulation of three-dimensional turbulent flow in the National Full-Scale Aerodynamic 
Complex (NFAC) Wind Tunnel Inlet 12 . 


TURBULENCE MODEL 

The prediction of turbulent flow in complex internal flow systems has now become 
realizable with the advent of supercomputers. One of the pacing items in this prediction 
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procedure is turbulence modeling. Turbulence modeling, especially at high speeds or in 
internal flows, is quite complex because of the paucity of appropriate experimental data 
and appropriate turbulence models, the two not being mutually exclusive. In internal 
flows, e.g., the boundary conditions at the inflow and outflow boundaries are crucial and 
not easily posed. In principle, some mean profiles and some turbulence profiles have to 
be imposed at the inflow boundaries in internal flows especially for elliptical flow fields. 
The sensitivity of turbulence models on the imposed inflow conditions in such situations 
is significant. 

Various approaches have been adopted to “close” the system of equations governing 
a turbulent flow-field. Algebraic models have been used with success for relatively simple 
flows, one or two-equation models for more complicated flow situations where large sepa- 
rated regions or severely strained flows exist. One popular model for internal flows is the 
two-equation, k - e model, which has been used with some success for massively separated 
and curved flows. However, different k - e models yield different results for the same test 
case. Various k - e models are different from one another through either wall function 
or low Reynolds number formulations. Away from the walls, all these models reduce to 
the so called “standard” k - e model. In spite of this fact, the predictions away from the 
walls differ for different models. One of the main reasons for this is that different models 
are validated with different experimental test data and since the experimental data also 
have some uncertainty about them, e.g., refer Sultanian et al 13 , the comparison of various 
turbulence models is not clearly made. However, direct numerical simulation results of 
turbulent flow in a channel, such as those of Kim et al 14 , can also be used to validate tur- 
bulence models, thus alleviating this experimental uncertainty problem. Another reason 
for the discrepancies in the model predictions, in the author’s experience, is that different 
inflow boundary conditions for the transport variables, k and e, are used in different models 
in integrating the k - e system. Since k and e are transport properties of the flow, their 
distribution throughout the flow field must be, to some exent, determined by the way they 
are prescribed at the inflow boundary in internal flows. 

In the present study, this particular problem is addressed and a low Reynolds number 
model is proposed which, 1) minimizes the sensitivity of the model to the manner in which 
the inflow boundary condition is prescribed, and 2) yields accurate predictions of the 
turbulence kinetic energy, k, and its peak near the wall. The present model is validated 
with reliable and detailed data 15 for a channel Reynolds number of 7,700 based on the 
channel height and also compared with the direct numerical simulation results of Kim et 
al 14 . The present model partly derives from some observations based on the results from 
direct simulation of near-wall turbulence, see Chapman and Kuhn 16 . 

In most of the low Reynolds number formulations, the dissipation rate at the wall is 
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set to zero and the effect of non-isotropy is incorporated through the governing equation 
for the kinetic energy, e.g., see Patel et al 17 . So the accuracy of near-wall predictions of 
these models is determined by the degree of approximation of the non-isotropic dissipation 
rate near the wall. For example, an inaccurate slope of the k - y curve near the wall and 
generally an under-prediction of the peak in the kinetic energy are a common feature in 
the predictions of these models. On inspection of the transport equations for k and c for 
a fully developed turbulent channel flow, it can be observed that the turbulence profiles 
are strongly dependent on the Reynolds number 18 . The inflow profiles must, therefore, 
appear as a set of curves with Reynolds number as a parameter. Therefore, in calculating 
a turbulent flow at a given Reynolds number, the inflow conditions must correspond to 
that particular Reynolds number. 

The present turbulence model has been used in the turbulence solver, KEM 5 , which 
solves the three-dimensional transport equations for the turbulence kinetic energy k and 
its dissipation rate e in generalized coordinates with the diffusion terms retained in all the 
three directions. It can be used in conjunction with any compressible or incompressible 
flow solver to compute the turbulent flows of interest. 


The governing equations for the transport quantities of turbulence, the turbulence 
kinetic energy, k and the dissipation rate, e, can be written in cartesian coordinates as the 
following 2x2 system 

d t D+(Fi-F Vi U=S 


with the turbulence kinetic energy and the homogeneous dissipation rate respectively given 

by 

k = -u l u l and e 




where the solution vector is 


D = 


pk 

pe 


the flux vectors are 


and the source term is 


Ft = 


pU t k 

and F„ = 

Pk^ii 

_pU t € 

* Re oo 



S = 


Rt r 


P - peRe oc - F 

C^P-Ctf^Re^-G 


The kinetic energy production term due to the mean shear is given by 

P = -pu l u j U lt] 
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where t/ t is the mean velocity and u t is the fluctuating velocity and < 5 t y is the Kronecker 
delta. The turbulent viscosity p t is given by 




CjxUlP}^ 

€ 




where C M is a constant and is a function that takes into account the low Reynolds 
number dependence of C M . 


In the present model, the function / M is defined as 

U = ^exp(-C 4 y + ) +exp(^-^) 

where C3 and C 4 are constants. The present prescription for is discussed later in the 
current section. 


Also, Hk = M + ^ and = p, + where cr* and <7 e , turbulent Prandtl numbers 
for the k and e transport processes respectively, are assumed to be close to unity. 

The function above takes into account the low Reynolds number dependence of the 
constant C2 and is defined 19 as follows: 

h = 1.0- ^|exp{ — (R r /6) 2 } 


The coefficient of viscosity p is normalized by p o©, velocities are normalized by U^, 
distances by a characteristic length L, and the density p by p^. This results in the reference 
Reynolds number definition as 

P P00U00L 

Moo 


The turbulent kinetic energy k is normalized by and the dissipation rate e is normalized 
by UHL. 


The equations for k and f used here are applicable for both high and low turbulent 
Reynolds number 



Values of various “constants” are chosen as follows: 

0 € 

Ok 

Ci 
C 2 

The values of C\ and C 2 are inferred 19 respectively from grid turbulence data and the 
requirement of consistency with the von Karman constant, k = 0.42. 

For the high Reynolds number form using wall functions as boundary conditions for 
k and e, the terms F and G are set to zero in the k- and e- equations respectively. To 
incorporate the low Reynolds number effects, the high Reynolds number form of the k- 
t system is modified in the present study as follows. The viscous diffusion of k and e 
is included in the governing equations. An appropriate term G is incorporated in the 
t equation, but the term F in the k equation is set to zero. This is different than the 
formulation used by Jones and Launder 20,21 and Chien 22 , where they set the term F equal 
to the value of the near-wall non-isotropic dissipation rate and then set the boundary 
condition, £ = 0 at the wall. The form of F is different in the two models 20-22 . However, 
since the differential equation for k involves the homogeneous dissipation rate and not 
just the isotropic dissipation rate, see, e.g., Hinze 23 , it is more appropriate to treat £ 
as the homogeneous dissipation rate and set the boundary condition on e at the wall 
corresponding to the non-isotropic part of the energy dissipation rate since its isotropic 
component reduces to zero at the wall. This boundary condition and the term G are found 
by using a Taylor series expansion procedure as follows. 

The three fluctuating velocity components close to the wall can be expressed as 

u = a x y + a 2 y 2 + a z y z + ... 

v = 6 2 t/ 2 + b z y z + ... 

and 

w = c x y + c 2 y 2 + c z y z + ... 

where a n ,b n ,c n ,n = 1,2,3 are functions of x, z and t, so that the turbulence kinetic 

energy near the wall is given by 

k — 2 [( a i 2 + c i 2 )y 2 + 2(aia 2 -f c x c 2 )y z + 0(y 4 )] (1) 
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and assuming streamwise and spanwise homogeneity near the wall, the dissipation 
rate near the wall is given by 


t = v 


(ai 2 + cj 2 ) + 4(axa 2 + cic 2 )y + 0(y 2 ) 


(2) 


Since the isotropic dissipation rate at the wall reduces to zero, the wall dissipation 
rate is entirely the non-isotropic component and it is given by 


f-w 



The zeroth order expression for the near-wall dissipation rate, e = which is in 

error by O(y), replaces the term F in an earlier formulation 22 . However, in another 
formulation 20 ’ 21 , the term F is expressed as 


which is a first order accurate expression for £. Either of the preceding two expressions 
for e can be used as a boundary condition for the dissipation rate at the wall. However, 
the former is more in the spirit of Dirichlet prescription and is more stable than the 
latter as experienced by the author 5,11 . An alternate first order accurate expression 
for e in terms of e w has been derived 4 and it is given by 


t P = \v — ^ - e w + 0(y 2 ) 

yp 2 

where P is a point sufficiently close to the wall. In their computations, Chapman and 
Kuhn 16 found that the resulting expression for the wall dissipation rate, e w , given by 
— ep varies only slightly upto y + « 1 and, therefore, they suggest that this 
boundary condition should be used when P is within y + « 1. The boundary condition, 
£ w = as used in Refs. 5 and 11, is also applied for y + ss 1. Combining the two, 

upto y + « 1, €p « Now, from the expression for k given above (Eqn. 1), it 

is clear that k varies as y 2 close to the wall. However, it is not clear how close to the 
wall this behavior is true. Since even in the most reliable experimental study (e.g., 
see Ref. 15), the fluctuating velocity components are not measured below y + « 1.5, it 
is reasonable to assume that k a y 2 in the region below y + « 1. As shown later, this 
prescription does yield the correct kinetic energy variation near the wall and a near- 
wall behavior of Reynolds stress, u'v', in qualitative agreement with that predicted 
by Coles 24 . 
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The form of the term G which is dissipation like in Ref. 22 balances the molecular 
diffusion term near the wall. In Refs. 20 and 21, the form of G is different and this term is 
added to match the kinetic energy with the experimental data at y + « 20. In Ref. 22, it is 
assumed that the length scale in the sublayer is proportional to the normal distance from 
the wall and that e <x y 2 ; therefore to be consistent with the Reichardt’s requirement 25 
that nt oc y 3 near the wall (see Refs. 14,16,24 for recent evidence of its validity for both 
streamwise homogeneous and inhomogeneous flows), a damping factor on the expression 
for /z t which goes to zero at the wall linearly with y has to be incorporated. 

In the present formulation, the term G in the e equation is such that the dissipation 
like term C2/2 ^R^oo is damped by a factor proportional to y 2 so that there is an effective 
balance between the two leading terms in the e equation near the wall , i.e., the molecular 
diffusion and the dissipation like terms. This is because e — » e w as y — ► 0; e w being of the 
order of 1. Accordingly, the term G is given by 

G = -C 2 fi?YRe 00 {exp(-C 4 y +2 )} 

The constants, C3 and C 4 are chosen so that the peak in the kinetic energy matches with 
the experimental value (Ref. 15). Their relative magnitudes are chosen so that the peak in 
the kinetic energy is located at the position given by the experiment (Ref. 15). It should 
be noted here that these constants are also tied to the type of the ramp function used 
for the length scale at the inflow boundary although the ramp function dependence of the 
present model is weaker than that of the other models 20-22 , as discussed in the section on 
Results and Discussion. 

Since the turbulent viscosity is given by \i t = and very near the wall, k <x y 2 

and (f - e w ) oc y, the function f ^ has been constructed so that it behaves as ^ very near 
the wall, (also see Ref. 16), which ensures the satisfaction of Reichardt’s requirement. 
Away from the wall, the function / M tends to 1. 

With this low Reynolds number form thus built into the k - e system, the governing 
equations can be solved subject to the boundary conditions applied at the wall. Before 
finite differencing the equations, they are re-cast in the generalized coordinate system. 
Using the generalized coordinate transformation, 

T = t 

Xj = Xj (*.-. t) 

where \j = (£,*?, c) an d x i = { x ,y,z), the k - e system transforms to 

d r D + d Xi (F t ~ F Vt ) - jS 
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where the solution vector, 


and the flux vectors, 



pk 

P* 


F 


1 

' pkQi' 

, -f 1 

ff s] J? 

MfcVx„V Xt d Xi k 

7 

peQi 

JR e<x 

. ^V Xl .Vx, d Xi e 


The Jacobian and the contravariant velocities are respectively given by 


J = 


d v 

-q 2 - and Qi = dtXi + Ujd Xj Xi 

x i 


Wall Boundary Conditions 

if >(. 

There are various expressions that can be derived for the near-wall and wall dissipation 
rates in terms of the kinetic energy. Some of these (e.g., Refs. 5,11,16,20-22) have been 
discussed above. Referring to Eq. 1 and Eq. 2, we can also obtain a zeroth and a first 
order expression for the near-wall dissipation rate and they are given by 


€ — 2i/ 


d(fc/y) 

dy 


+ 0(y 2 ) 


and 


1 dk , 
t = V ydy +0(y) 


These expressions can also be used as boundary conditions on e. However, the boundary 
condition, 


£ W 



or 


iw - 4u— - £ + 0(y 2 ) 


are the most attractive since they do not involve computing any gradients and therefore 
are more stable. 
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Inflow Boundary Conditions 


Since the kinetic energy is assumed to vary parabolically upto y + « 1 , therefore, let 
k = ay 2 in that region. The constant a is calculated as a = where kp is taken from 
experiment at a point P distant yp from the wall corresponding to y + « 1. In this region, 
the length scale is given by / a y 3 / 2 as opposed to / = /cy in the buffer or overlap layer. 
This follows from the mixing length hypothesis, 


Pt = pi 2 


dU 

dy 


and the fact that p t <x y 3 very near the wall. The constant of proportionality in the 
expression for the length scale is calculated by matching the two expressions, / a y 3 / 2 and 
l — Ky , at y+ « 1. Therefore, if / = /?yp 3/2 , then ft = ^=, where y P corresponds to 

y+ ss 1. 

, 2 ^ 

On comparing p t = p/ 2 and we have an expression for the dissi- 

pation rate very near the wall as 


* = 


k 2 



The value of |^| is the slope of the experimental curve for U vs. y, which is a straight 
line throughout the sublayer. This curve was extracted from the information in Figs. 3 
and 4 of Ref. 15 and the values therein for u T and U at the centerline. Away from the 
wall, the dissipation rate is given by 


C„ 3 / 4 fc 3 / 2 


The preceding specification on / and e is used to construct the inflow boundary con- 
dition on e which is then fixed during the course of computations. This obviates the need 
for specifying a ramp function for / in the near-wall region to construct an e profile unlike 
in the other models. This is important since otherwise the predictions (as in other models) 
suffer progressively as the grid near the wall is refined. In the course of the present com- 
putations, it was observed that for a given model it is the cut-off value of the length scale 
ramp function that determines the value of k away from the walls and at the center-line. 
As the slope of the / - y curve below the cut-off point is changed, this cut-off value has also 
to be changed to match the kinetic energy away from the wall, i.e., if the slope of the / - 
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y curve is increased, the cut-off value is decreased. In the present model, the predictions 
are seen to be less sensitive to the slope of ramp function since it is used above y + « 1. 
However, in the other models, the smaller the slope of the / - y curve is prescribed, the 
larger should this cut-off value be, and the peak in the kinetic energy is not matched as 
accurately as before. 


CURVATURE CORRECTION 

2 

In the curvature model 1 , the decay term, C2/2 in the e equation is multiplied 

by (1.0 — C c Rit), where the rotational Richardson number, Rit = V " c r ° sa ^§7^ and 
the coefficient, C c , is a model parameter. In the study of turbulent flow over spinning 
cylinders 1 , a constant value in the range, 0.1 - 0.2 (depending on the application), was 
used. In the present study, this parameter is made to vary both along the angular and 
crossflow directions. This form is preliminary and may need some modification when 
extensive experimental data are available for validation of the model. For the present 
study, C c is set to 0.1 from 0° to 45°, to 0.2 from 45° to 135° and back to 0.1 from 
135° to 180°. A dependence of this parameter on the crossflow direction is also indicated. 
Various linear variations from the convex to the concave walls were tried; a monotonic 
dependence of the solution on the magnitude of this variation was observed. Fine-tuning 
of this parameter will be achieved after a detailed comparison with the experimental data 
is made. 


NUMERICAL METHOD 


The k - c system in generalized coordinates is integrated using the implicit noniterative 
approximate factorization algorithm of Beam and Warming 10 . The essential details of the 
scheme are given below. Before finite differencing the governing equations, the flux vector 
F Vi is redefined as 


1 

Pk{pk) „ 

1 

pkpk{p ) ,, 

pRt 0 o 

./*<(/*) H . 

8 

0$ 

1 

. PtPe{p) H . 


so that the second term is lumped explicitly with the source term, S. Writing the finite 
difference expression for the governing equation in generalized coordinates, we have 


D n+1 

J»+i 


+ At 


- r*) 


n+ 1 




where the time differencing is first-order (Euler) implicit, and the spatial derivatives are 
approximated by the central-differencing operator 6 Xi . 
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Linearizing the flux vector F t locally in time, this equation can be written as 
{/ + J n+l A t f Xl (m* - } (D n+1 - D") = (D n + S At) - ID n - 

A^ n+, ft, l(F. n ~ - p 6 Xi D-y+^ (< X1 P)C”|} 


where the Jacobian matrix of the flux vector F, is given by 


M i = (A,B,C)=d- 5 F i n = ± 


Qi 0 
0 Q t 


with 


a = 


r 

J n +'Re 0 o 


Vx.Vx, 


r = 


Hk o 
0 


The right hand side of this equation is calculated explicitly, and the left-hand side can be 
factored approximately. Therefore, we have 


/ + ArJ n+1 ^ ^4 n 
1+ ArJ n+1 6s ( C n 


+ c.V* A* I + A Tj n+l Sr, (V - - p 6^j + e t V„A r 

- p 6^ + 6<V f Aj ( D n+1 - D n ) = RHS - t e [(Vx.Ax*) 2 ] D n 


where V and A are first order forward- and backward- difference operators, RHS is the 
right hand side of equation (4), e, is the implicit second order smoothing parameter, and 
€ e is the explicit fourth order smoothing parameter. The smoothing terms are added to 
damp the high frequency oscillations in the solution which arise out of the central-difference 
approximation of the spatial derivatives. 


DISCUSSION OF RESULTS 

The experimental data used for validating the low Reynolds number turbulence model 
is that of Kreplin and Eckelmann 15 . The Reynolds number based on the channel height 
is 7,700. The validation data corresponds to fully developed conditions. Therefore, the 
computations were also carried out at the fully developed mean flow condition. Fully 
developed mean velocity profile was imposed throughout the computational domain from 
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the experiment and only the transport equations for k and e were integrated. This simpli- 
fies the problem considerably in terms of validation of the turbulence model because the 
algorithm dependence of the mean velocity field solver is entirely eliminated. 

Calculations were carried out for a coarse and a fine grid; in the former, the first point 
from the wall was situated at y + « 2.0 and in the latter, it was located at y + ss 0.1 from 
the wall. This was necessary so as to check whether the model was uniformly valid as 
y + — *• 0. The present predictions are compared with those of Jones and Launder 20 ’ 21 and 
Chien 22 . The results corresponding to these models were obtained using the KEM code 5 
in the same fashion as the present model. 

Given the turbulence model, a distribution for the length scale was specified for the 
t prescription at the inflow boundary that would yield a k profile in best agreement with 
the experiment, especially the peak in the k profile. A ramp function / = k y, y < 0.35; / = 
constant, y > 0.35 was used for the earlier two models 20-22 as well as the present model. 
This resulted in the best prediction of the peak in kinetic energy in all the three models. 
For simplicity, this length scale distribution will be called distribution A. 

An alternative distribution for the length scale, referred to as distribution B, was 
prescribed as follows. A ramp function / = 2.5 y,y < 0.0136;/ = constant, y > 0.0136 
was used for the model 20-21 *, for the model 22 , a ramp function / = 2.5 y,y < 0.033;/ = 
constant, y > 0.033 was used. For the present model, a ramp function, / = 2.5y, y < 
0.016;/ = constant, y > 0.016 was used. 

Results corresponding to the length scale distribution A are discussed below. The 
results corresponding to distribution B will be discussed next via a comparison of kinetic 
energy predictions. 

In Fig. 1, the results are shown corresponding to the coarse grid in terms of the kinetic 
energy variation from the wall to the center-line. The kinetic energy is normalized with 
respect to u T 2 ; the friction velocity, u T , is fixed from the experiment. The experimental 
peak in k/u T 2 is matched by the present model with no visible discrepancy whereas the 
two earlier models do not predict the peak very accurately. The slope of the k - y curve 
near the wall as predicted by the Jones and Launder model 20,21 matches that due to the 
present model because the near-wall non-isotropic dissipation rate has been approximated 
by a first order expression and not a zeroth order expression as in the Chien’s model 22 , 
are in good agreement with each other. The dissipation rate profiles are shown in Fig. 2. 
The e profile exhibits a linear behavior near the wall with the present model and e attains 
a non-zero value at the wall while it goes to zero at the wall in the other two models. 
However, implicit in the t profiles corresponding to the latter two models is the absence 
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of the near-wall non-isotropic dissipation rate. 

Results corresponding to the fine grid are shown in Fig. 3 and Fig. 4, and they 
show a comparison among the experiment 15 , the present model and Jones and Launder 
model 20 ’ 21 . As shown in Fig. 3, the experimental peak in kinetic energy is accurately 
matched by the present model, whereas the latter model underpredicts it appreciably. 
This shows that the present model is uniformly valid as y + — ► 0. The corresponding e 
profiles are shown in Fig. 4. 

A comparison among the experimental data, direct simulation results of Kim et al 14 
and the present predictions is shown in Fig. 5. The kinetic energy is plotted against the 
wall variable, y + — J ^ L . The predictions shown correspond to the coarse grid. Since the 
present model is validated with the experimental data and not the direct simulation results, 
the agreement between the present predictions and the results from the direct simulation is 
not as good. For an added comparison, the near- wall prediction of the normalized Reynolds 
stress, - , has been compared with that corresponding to Coles’ deterministic model 

(see Ref. 24) which gives the Reynolds stress as 

uV - 0.00103{^} 3 + 0(y 4 ) 

The Reynolds stress is plotted versus y + in Fig. 6. The behavior of the present predictions 
is in reasonable agreement with that of Ref. 24. If the model was validated with the 
direct simulation results of Kim et al 2 , it is reasonable to expect that this Reynolds stress 
comparison would be better because the level of turbulence near the wall from the direct 
simulation is not as high as in the experiment. 

The behavior of Reynolds stress is strongly dependent on the present / M prescription 
and this prescription is based on the model validation with respect to the experimental 
data 15 . It is also interesting to note that in the experiment 15 , the viscous and the turbulent 
stresses become equal at y + ss 13, and that in the present computations, the viscous stress 
equals the Reynolds stress, - u'v' , at y + « 11 in the coarse grid case and at y + « 12 in 
the fine grid case. 

The kinetic energy variation from the wall to the center-line corresponding to length 
scale distribution B is shown in Fig. 7. All the three models predict the peak very 
accurately. The experimental data at the center-line are matched accurately only by the 
present model. These predictions correspond to the coarse grid case. 

Computations for turbulent flow in the curved duct have been carried out in two 
stages. First, an approximately fully developed two-dimensional turbulent channel flow 
profile at the outflow boundary was obtained by running a straight channel case at a 
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Reynolds number of 89,160 based on the channel height. The inflow boundary conditions 
for this case were taken from experiment 26 ; the inflow profiles for the mean velocity and the 
turbulence kinetic energy were thus prescribed. In the second stage, the inflow boundary 
condition for the curved duct was fixed to be the same as the outflow boundary condition 
from the first stage of the computations for a straight channel. 

The grid for the curved duct was chosen to have 99 points across its width (in the 
crossflow direction) and 141 points in the streamwise direction. The outflow boundary 
was placed about 25 duct heights downstream of the inflow boundary. The grid for the 
U-bend is shown in Fig. 8. The radius of curvature of the curved duct is the same as 
for the Space Shuttle Main Engine (SSME) Turn-Around Duct (TAD), and is given by 
ri/H = 0.5, where r t is the radius of the inner (convex) wall and h is the duct height. 

The streamfunction contours in the U-bend are shown in Fig. 9. The separation 
zone on the convex wall extends from approximately the 120° location to over a duct 
height downstream of the 180° location. In Fig. 10, mean streamwise velocity is plotted 
for various angular locations from the convex to the concave walls. The flow is seen to 
accelerate upto slightly upstream of 90° location near the convex wall and then it starts 
to decelerate downstream. It eventually separates at about 120° location and reattaches 
over a duct height downstream of 180° location. This is in good agreement with the 
experimental results 27 obtained at Colorado State University. 

The distribution of the turbulence kinetic energy in the U-bend is shown in Fig. 11 
through the k contours. High levels of kinetic energy are seen near the concave wall 
especially across the separation bubble on the convex wall. Higher levels of kinetic energy 
are seen all around the bend from 0° to 180° on the concave wall as compared with those 
on the convex wall. The location of sharp gradients in the kinetic energy profiles across 
the duct is seen to be confined close to the wall on both sides as seen in Fig. 12. These 
peaks in the kinetic energy distribution represent high levels of turbulence production due 
to the presence of counter-rotating vortices close to the wall as seen in the case of a straight 
channel. Beyond this peak, away from the walls, higher levels of kinetic energy are present 
near the concave wall. This is due to the effect of curvature on turbulence. Since this 
study is only two-dimensional, no Goertler vortices are taken into account. Therefore, 
higher levels of energy near concave wall do not reflect any effect of Goertler vortices, but 
just the effect of rotational Richardson number introduced through the decay term in the 
e equation. Generally, the level of kinetic energy near the concave wall is seen to be two to 
three times larger than that on the convex wall. Fig. 13 shows the variation of turbulence 
viscosity from convex to concave walls for various angular locations along the bend. Higher 
levels of turbulence viscosity are present near the concave wall. In general, the behavior 
of turbulence viscosity reflects that of the turbulence kinetic energy. Turbulence shear 
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stress profiles are shown in Fig. 14. The change in magnitude of shear stress from convex 
to concave walls is about an order of magnitude at, e.g., 90° and 128° locations. This 
order of magnitude variation of shear stress was also reported by Rocketdyne 28 in their 
experimental study of turbulent flow in an axisymmetric curved duct. 

Skin friction on the two walls is shown in Fig. 15. Skin friction on the inner wall, C/ ; , 
shows the extent of the recirculation bubble on the convex wall. Some data points from 
Monson’s experiment 29 conducted at NASA Ames Research Center at a Reynolds number 
of 10 6 based on the duct height are also plotted in Fig. 15. Since a sufficient number 
of data points are not available, the comparison with experiment is still preliminary. It 
is reasonable to expect that the peak in experimental C/ ; should lie between S/(H/2) = 
0.56 and 2.16, see Fig. 15. However, this peak should be lower than the peak predicted 
by the computations since the Reynolds number corresponding to the computations is 
approximately 90,000, whereas that for the experiment is 10 6 . That the peak in Cf i is 
upstream of the start of the bend is quite interesting. It would appear from the mean 
velocity profiles that the skin friction on the convex wall should peak out around 9 = 52° 
due to a sharp velocity gradient there, but the turbulence viscosity at the corresponding 
location is lower than that at those locations upstream of the bend and hence the shift of 
the peak in C j t upstream of the 52° location. This is a useful piece of information for the 
experimentalist who can thus look for and measure the quantities of interest accurately and 
efficiently. Although the flow conditions and the geometry are not the same, the behavior 
of C f ti , skin friction on the outer wall, predicted by the present computations is similar to 
that shown in Fig. 10 in Ref. 30. The skin friction on the outer wall does not increase 
monotonically, but there is a slight jitter about some mean monotonic variation as shown 
in Fig. 15 and in Fig. 10 in Ref. 30. The asymptotic values of C/ t and C/„ do eventually 
approach the fully developed turbulent channel values corresponding to Re a = 90,000. 

Finally, a comparison is shown between the present two-dimensional predictions and 
the axisymmetric experimental data of Sharma and Ostermier 28 in Fig. 16. Although the 
present computations are carried out in a two-dimensional formulation and not axisym- 
metric, the trend in the velocity profiles predicted by the computations is in qualitative 
agreement with that given by the axisymmetric experiment. 

CONCLUDING REMARKS 

A low Reynolds number form of the k - e turbulence model has been presented. The 
calculations were carried out by prescribing the fully developed mean velocity profile from 
experiment throughout the solution domain and by only solving for the turbulence quan- 
tities. This procedure eliminates the uncertainties that creep into the model predictions 
through the algorithm dependence of the mean flow solver. The model has been validated 
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with experimental data and has been evaluated with the results from the direct simulation 
of turbulence in a channel. The model predictions are also compared with those of earlier 
models. 

The model predictions of the near-wall turbulence kinetic energy are in very good 
agreement with the experimental data. The experimental peak in the kinetic energy is 
matched exactly. The predictions are in reasonable agreement with the direct simulation 
results. The model is less sensitive to the specification of the inflow boundary condition on 
t than other earlier models. However, there is a definite need for the upstream experimental 
data at the inflow boundary in terms of the turbulent kinetic energy appropriate to the 
particular Reynolds number of interest. Any approximations made in specifying the inflow 
boundary conditions will correspondingly influence the solution downstream especially 
when using a low Reynolds number turbulence model. 

The function, / M , is so prescribed that the model is validated with the experimental 
data. Therefore, the near-wall behavior of Reynolds stress as predicted by a deterministic 
model is not accurately matched by that predicted by the present model. The peak in 
the kinetic energy is predicted accurately for both the coarse and fine grids. By refining 
the grid, little change in the predictions of kinetic energy was seen to occur and therefore 
not only is the solution grid independent but the turbulence model is uniformly valid as 
y + ->0. 

Alternate boundary conditions for c have been derived although not used in the present 
calculations. 

The present study also deals with the simulation of turbulent flow in a U-bend using 
a low Reynolds number k - c turbulence model. The turbulence model includes an existing 
curvature correction which has been slightly modified in the present study. Detailed pre- 
dictions on mean streamwise velocity, turbulence kinetic energy, skin friction, turbulent 
viscosity, etc., are made. The predicted results are in agreement with the physics of the 
problem. Some quantitative comparison is made with the available experimental data. 
An approach to refine the existing curvature correction model used here is presented that 
takes into account the dependence of the curvature parameter on the space variables. 

The predictions reveal features of the flow field that can serve as a guide in fine-tuning 
experimental investigations of this class of flow fields. The present turbulent flow solution 
capability can be used in many engineering design situations such as in the Space Shuttle 
Main Engine (SSME) Turn-Around Duct (TAD) turbulent flow processes. 
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Figure 3. Distribution of turbulence kinetic energy from the wall to the 
center-line for a fine grid 
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Distribution of dissipation rate from the wall to the center-line 
for a fine grid 
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Figure 6. Near-vvall variation of Reynolds stress with the wall variable, y + , 
for a fine grid 
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Figure 7. Distribution of turbulence kinetic energy from the wall to the 
center-line for a coarse grid and length scale distribution B 






Figure 8. Grid for a Two-Dimensional U-Bend (180° Curved Duct) 
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Figure 12. Cross-Flow Variation of Turbulence Kinetic Energy at Various 
Angular Locations 
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Figure 13. Cross-Flow Variation of Turbulence Viscosity at Various Angular 
Locations 
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Figure 16 . Comparison Between Predictions and Experimental Data for Mean 
Streamwise Velocity at Various Sections Around the Bend 
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Figure A-l. Convergence History of the KEM Code Corresponding to Wall 
Function and Wall Boundary Condition Approaches 
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APPENDIX 


User Manual for KEM Code 

The code, KEM, has been developed in a fashion parallel to that of the AF algorithm 
codes developed at NASA Ames Research Center; consequently, it can easily be interfaced 
with any of those codes for simulating complex turbulent flows. Many of the integer 
parameters and counters used in the incompressible flow code, INS3D, are similarly used 
in KEM. There are some additional parameters used in KEM that are given default values 
in the initialization routine, INITKE, but are updated through the NAMELIST statement. 
The turbulence model constants are initialized in the routine INITKE. 

A short description of each individual routine follows. 

ENERGY 


This routine is the main routine in the code, KEM, which is interfaced with the code 
INS3D, so that all the relevant parameters and variables are transferred back and forth 
from one to the other. The initialization routine, INITKE, is called here and the array 
VNUT(I) is rescaled as it is transferred from the INS3D. After VNUT(I) is updated at the 
end of one complete sweep in KEM, it is rescaled back to be compatible with INS3D. The 
routine STEPKE is called here, which is the driver routine that performs all operations 
necessary to update the arrays QT(1,I) and QT(2,I), pk and pe respectively, in the interior 
of the solution domain from one time step to the next. Immediately after calling STEPKE, 
ENERGY calls BCKE. the routine that updates the variables at the boundaries. Following 
this, the array VNUT(I) is calculated with the new values of pk and pe, and the the new 
values of arrays, QT(1,I), QT(2,I) and VNUT(I), are printed through the routine OUTPT 
and stored, as and when desired. 

INITKE 


This is the subroutine to input data and initialize various parameters in the turbulence 
model. It also initializes the arrays QT(l,I), QT(2,I) and VNUT(I). It calls CHFTKE, 
which fits the experimental data, smoothes the data by calling SMTH, and thereby provides 
the inflow boundary conditions on pk and pe. The routine INITKE then initializes the 
variables throughout the solution domain. 
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STEPKE 


This driver routine performs functions in the same manner as the routine STEP in 
INS3D. It calls RHSKE, which calculates the right hand side of the pk — pt system of 
equations and thereby updates the arrays ST(1,I) and ST(2,I). Then it calls the routines, 
FLTRKT and BLTR2, which accomplishes the matrix inversion in the K-direction of the 
approximate factorization algorithm and thereby the array ST is updated. The routines 
PBLTR2 and BLTR2 accomplish the 2x2 block tri-diagonal inversion in a given direction 
respectively depending on whether the given direction is periodic or non-periodic. Con- 
ventionally, in both INS3D and KEM, it is the K-direction that is allowed to be either 
periodic or non-periodic. 

Similarly, the inversions in the J- and L-directions are accomplished respectively by 
first calling FLTRJT and BLTR2 and thereby updating the array ST and then by calling 
FLTRLT and BLTR2 and then updating the arrays ST and QT. 

RHSKE 

This routine calculates the right hand side of the implicit pk - pt system. The routines 
METRIC, FLUXV and VSRHST are called for each of the three J-, K- and L-directions, 
and respectively return the metric quantities, the inviscid flux vector and the viscous right 
hand side terms appropriate to the given direction. Having thus partially updated the 
array ST, it calls the routine SMOOKE which calculates the explicit smoothing terms and 
lumps them into the array ST. Finally, the routine SOURCE is called which calculates the 
source terms in the pk - pt system and thereby updates the array ST. 

FLUXV 


This routine calculates the inviscid flux vector of pkU and ptU for the J-direction, 
pkV and ptV for the K-direction and pkW and ptW for the L-direction; U, V and W are 
the corresponding contravariant velocities. 

V SRHST 

This routine calculates the viscous terms of the right hand side completely locally for 
all of the three co-ordinate directions. The density terms corresponding to the compressible 
flow fromulation are also calculated here. The array ST is thereby updated in this routine. 
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SMOOKE 


This routine adds the explicit smoothing terms to the right hand side of the pk — pe 
system corresponding to each of the three directions. All the calculations are done locally. 
The array ST is returned to the calling routine RHSKE. 

SOURCE 


This routine is called by RHSKE at the end. All the calculations are done locally. 
Most of the effort is expended in computing the appropriate metric terms in the three 
directions. However, if the metric quantities were to be calculated initially and stored for 
future use, the additional expense of having to calculate the metric terms in SOURCE 
would be avoided. The array ST is updated and returned to RHSKE. 

FLTRJT, FLTRKT, FLTRLT 

These routines fill the 2x2 blocks AT, BT and CT of the tri-diagonal matrix in each 
of the three co-ordinate directions by calling METRIC and AMATX corresponding to the 
appropriate direction. The implicit smoothing terms are also included in the arrays AT, 
BT and CT. 

AMATX 

The Jacobian matrix of the flux vector is computed here corresponding to a given 
direction and the array AT is returned to the appropriate calling routine. 

BCKE 


This routine imposes the appropriate boundary conditions. The wall function formula- 
tion is used if the logical variable, IWAL, is defined to be .TRUE, through the NAMELIST 
statement; otherwise, the low Reynolds number formulation is used. All the calculations 
are done locally for the latter case; in the former case, the routine LOGLAW is called 
which iteratively calculates the friction velocity, VTAU. 

LOGLAW 


The friction velocity, VTAU is calculated here iteratively and the array VTAU is 
returned to BCKE after having converged to within a specified tolerance. 
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BLTR2 (PBLTR2) 


This routine accomplishes the inversion of tridiagonal system of 2x2 block matrices 
corresponding to each of the three coordinate directions, (for periodic K-direction PBLTR2 
is called), and returns the solution vector F to the calling routine STEPKE. The routine 
LUDEC2 is called here. 

LUDEC2 


This routine computes the L - U decomposition elements. Block inversions in the 
routine BLTR2 (PBLTR2) use non-pivoted L-U decomposition. 

OUTPT 


The indices I, J, K and L, and arrays QT(1,I), QT(2,I), VNUT(I), Q(N,I), X(I), Y(I), 
Z(I) and VNUT(I) are printed for a given co-ordinate direction. This routine is called by 
ENERGY when NTIME equals NPRNT. 


There are some default parameters that are set in the INITKE routine that can be 
conveniently over-ridden through a NAMELIST statement in INITKE and then assigning 
the new values to them at the end of the source code. These parameters are listed below. 

NAMELIST/D ATAKE/'DISKEN, IWAL, JOFF, KOFF, LOFF, 

SMUK. SMUIMK 

The first two parameters DISKEN and IWAL are defined to be logical. DISKEN 
is .TRUE, for a restart run, and in that case, a stored solution for the arrays QT and 
VNUT is read in ENERGY. If IWAL is .TRUE., then the wall function formulation is 
used in BCKE. The parameters JOFF, KOFF, LOFF, are set to zero if IWAL is .FALSE.; 
otherwise, appropriate values are given them, depending on how far away from the walls 
the wall functions are prescribed. The wall function formulation is tested out for KOFF 
.NE. 0; the other two directions have not been tested yet. The parameters SMUK and 
SMUIMK respectively are the explicit and implicit smoothing parameters. 


Sample Problems 

Amongst other applications, the k - e turbulence model was used to compute the 
turbulent flow in a channel and a duct; the KEM code contains these two options. For 
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simplicity, the side-wall boundary conditions for the duct are taken corresponding to those 
of a symmetry plane. The computational mesh used is a Cartesian grid system which is 
stretched in streamwise and crossflow directions. The k - e solver, KEM, is made to lag the 
INS3D flow solver by one time step. To test the convergence property of the turbulent flow 
solver, a known, fully developed turbulent channel velocity profile was imposed and the k 
- e solver was iterated to produce a converged solution of the k - c system by switching 
off the INS3D. The convergence history is shown in Fig. A-l, where the root-mean- 
square values of the change in the solution vector from one time step to another is plotted 
against the normalized time r. These root-mean-square values are denoted by RMSD. A 
comparison is made between the wall boundary condition and the wall function approach. 
The wall function approach exhibits superior convergence characteristics. However, both 
wall boundary condition and wall function approaches yield a very fast convergence on the 
k - e solution as shown by the RMSD plots. 
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